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Abstract 



Iterative Rational Krylov Algorithm (IRKA) of [8] is an interpolatory model reduc- 
tion approach to optimal 'H2 approximation problem. Even though the method has 
^^ , been illustrated to show rapid convergence in various examples, a proof of convergence 

^/■^ ' has not been provided yet. In this note, we show that in the case of state-space sym- 

^^ . metric systems, IRKA is a locally convergent fixed point iteration to a local minimum 

of the underlying ?^2 approximation problem. 



1 Introduction 



k> ! Consider a single-input-single-output (SISO) linear dynamical system in state-space form: 

x{t) = Ax{t) + bu{t), y{t) = c^x{t), (1) 

where A G M"^", and b,ce M". In ([T]), x{t) G M", u{t) G M, y{t) G M, are, respectively, the 
states, input, and output of the dynamical system. The transfer function of the underlying 
system is H{s) = c^{sl — A)~^b. H{s) will be used to denote both the system and its 
transfer function. 

Dynamical systems of the form ([T]) with large state-space dimension n appear in many 



applications; see, e.g., [l[ and [10[. Simulations in such large-scale settings make enormous 
demands on computational resources. The goal of model reduction is to construct a surrogate 
system 

Xr{t) = ArXr{t) + brU{t) , yj.{t) = cJXr{t), (2) 



of much smaller dimension r '^ n, with Aj. G R''^'" and br, c,. G K*" such that yr{t) ap- 
proximates y{t) well in a certain norm. Similar to H[s), the transfer function Hr{s) of the 
reduced-model ([2]) is given by Hr{s) = cJ{sIr — Ar.)~^br. We consider reduced-order models, 
Hr{s), that are obtained via projection. That is, we choose full rank matrices Vr, Wr G R"^'' 
such that W^Vr is invertible and define the reduced-order state-space realization with (j2]) 
and 

Ar = {W,^Vrr^W^AVr, br = {W^Vry^W^b, Cr = V^ C. (3) 

Within this "projection framework," selection of Wr and Vr completely determines the 
reduced system - indeed, it is sufficient to specify only the ranges of Wr and Vr in order to 
determine Hr{s). Of particular utility for us is a result by Grimme |6|, that gives conditions 
on Wr and Vr so that the associated reduced-order system, Hr{s), is a rational Hermite 
interpolant to the original system, H{s). 

Theorem 1.1 (Grimme [6|). Given H(s) = c'"(sI — A)^^b, andr distinct points Si, . . . ,Sr G 
C, let 



Vr = [{sil -A)-^b... {srI - A 
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(4) 



Define the reduced-order model Hr{s) 
Hermite interpolant to H at Si, . . . ,Sr 



zKsl 
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Hr{s,) 



and 



Ar) ^br as in ^ Then Hr is a rational 
H'{si) = Hl{si) for z = l,...,r. (5) 



Rational interpolation within this "projection framework" was first proposed by Skelton et 



method of Ruhe 14 



al. P^ , [20[ , [21| . Later in p], Grimme established the connection with the rational Krylov 



Significantly, Theorem 11.11 gives an explicit method for computing a reduced-order system 
that is a Hermite interpolant of the orginal system for nearly any set of distinct points, 
{si, . . . , Sr}, yet it is not apparent how one should choose these interpolation points in order 
to assure a high-fidelity reduced-order model in the end. Indeed, the lack of such a strategy 
had been a major drawback for interpolatory model reduction until recently, when an effective 
strategy for selecting interpolation points was proposed in [8[ yielding reduced-order models 
that solve 
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where the H2 system norm is defined in the usual way: 
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The optimization problem ([6]) has been studied extensively, see, for example, 13|, |l9|, |8|, [l5|, 
y, [iTI, |7|, y, y, |22| and references therein. ([6]) is a nonconvex optimization problem and finding 



global ininiinizers will be infeasible, typically. Hence, the usual interpretation of ([6]) involves 
finding local minimizers and a common approach to accomplish this is to construct reduced- 
order models satisfying first-order necessary optimalitv conditions. This may be posed either 



in terms of solutions to Lyapunov equations (e.g., 19|, ll5|, l22|]) or in terms of interpolation 
(e.g., B a 03): 



Theorem 1.2. (I la . \BJ) Given H{s), let Hr{s) be a solution to ^ with simple poles 
Ai,...,Ar. Then 

H{-X.i) = Hri-X.i) and H'{-X.i) = Hl{-Xi) for z = l,...,r. (8) 

That is, any 'H2-optimal reduced order model of order r with simple poles will be a Hermite 
interpolant to H{s) at the reflected image of the reduced poles through the origin. 

Although this result might appear to reduce the problem of ?^2-optimal model approxima- 
tion to a straightforward application of Theorem 11.11 to calculate a Hermite interpolant on 
the set of refiected poles, {— Ai, . . . , — A.^}, these pole locations will not be known a priori. 
Nonetheless, these pole locations can be determined efficiently with the Iterative Rational 
Krylov Algorithm (IRKA) of Gugercin et al. [8|. Starting from an arbitrary initial selection 
of interpolation points, IRKA iteratively corrects the interpolation points until ([H]) is satisfied. 
A brief sketch of IRKA is given below. 



Algorithm IRKA. Iterative Rational Krylov Algorithm 

Given a full-order H{s), a reduction order r, and convergence tolerance tol. 

1. Make an initial selection of r distinct interpolation points, {sj}^, that is 
closed under complex conjugation. 

2. Construct Vr and Wr as in (jlj). 

3. while (relative change in {si} > tol) 

a.) Ar = iW,^Vr)-^W^^AVr 

b.) Solve r X r eigenvalue problem AjU = Au and assign Sj < Xi{Ar) for 

z = l,...,r. 

c.) Update Vr and Wr with new Sj's using (jl]). 

4. Ar = {W^Vr)~^W^AVr, K = {W^Vr)-^W^b, and Cr = V/c. 



IRKA has been remarkably successful in producing high fidelity reduced-order approximations 
and has been successfully applied to finding 'H2-optimal reduced models for systems of high 
order, n > 160,000, see (9[. For details on IRKA, see |8[. 



Notwithstanding typically observed rapid convergence of the IRKA iteration to interpolation 
points that generally yield high quality reduced models, no convergence theory for IRKA has 
yet been established. Evidently from the description above, IRKA may be viewed as a fixed 
point iteration with fixed points coinciding with the stationary points of the I-L2 minimization 
problem. Saddle points and local maxima of the I-L2 minimization problem are known to 



be repellent [ll[. However, despite effective performance in practice, it has not yet been 



established that local minima are attractive fixed points. 

In this paper, we give a proof of this for the special case of state-space-symmetric systems 
and establish the convergence of IRKA for this class of systems. 



2 State-Space-Symmetric Systems 

Definition 1. H{s)=c^{sl — A)^^b is state- space- symmetric (SSS) if A=A^ and c = b. 

SSS systems appear in many important applications such as in the analysis of RC circuits 
and in inverse problems involving 3D Maxwell's equations [J]. 

A closely related class of systems is the class of zero-interlacing-pole (ZIP) systems. 

n-l 

Definition 2. A system H{s) = K^ is a strictly proper ZIP system provided that 

n (^ - A,) 

> Ai > 2:1 > A2 > -22 > A3 > ■ • • > Zn-l > Xn- 

The following relation serves to characterize ZIP systems. 

Proposition 2.1. tl a] H{s) is a strictly proper ZIP system if and only if H{s) can be written 

as H{s) = 2_] — —r with Aj < 0, bi > 0, and Aj 7^ Xj for all i 7^ j. 
^ s — Xi 

The next result clarifies the relationship between SSS and ZIP systems. 

Lemma 2.1. ll^ J Let H(s) be SSS. Then H(s) is minimal if and only if the poles of H{s) 
are distinct. Moreover, every SSS system has a SSS minimal realization with distinct poles, 
and is therefore a strictly proper ZIP system. 

It can easily be verified from the implementation of IRKA given above, that for SSS systems, 
the relationship Vr=Wr is maintained throughout the iteration, and the final reduced-order 
model at Step 4 of IRKA can be obtained by 

Ar = QjAQr br = Cr = Q^b, (9) 

where Qr is an orthonormal basis for V^; the reduced system resulting from IRKA is also SSS. 



3 The Main Result 

Theorem 3.1. Let IRKA be applied to a minimal SSS system, H{s). Then every fixed point 
of IRKA which is a local minimizer is locally attractive. In other words, IRKA is a locally 
convergent fixed point iteration to a local minimizer of the ^2 optimization problem. 

To proceed with the proof of Theorem 13.11 we need four intermediate lemmas. The first 
lemma provides insight into the structure of the zeros of the error system resulting from 
reducing a SSS system. 

Lemma 3.1. Let H{s) be a SSS system of order n. If Hr{s) is a ZIP system that interpolates 
II{s) at 2r points Si, S2, . . . , S2r, not necessarily distinct, in (0,oo), then all the remaining 
zeros of the error system lie in (— oo,0). 

Proof. By Lemma 12.11 we may assume that II{s) is a strictly proper ZIP systems. Since 
H{s) is a strictly proper ZIP system, its poles are simple and all its residues are positive. Let 
Aj < 0, 0i > 0, for z = 1, . . . , n be the poles and residues of -ff(s), respectively. Now let 

2r n—T—1 n r 

Ris) = l[{s-s,), P{s)= l[{s + z,), Q{s) = l[{s-X,), Q{s) = l[{s-~X,), 

i=l i=l i=l i=l 

where A,, s,, and Zi are, respectively, the poles of Hr{s), the interpolation points, and the 
remaining zeros of the error system. Then for some constant K, H{s) — Hj.{s) = K J)'^'-^'^' . 

First suppose that {Xi}^^i fl {Afc}^^^ = 0. Then for each Xj, j = 1, . . . ,n, 

ResiHis) - Hris); A,) = K_ZM^i^z) = </*,> 0. (10) 

n (A, - A,)Q(A,) 

i=l 

Thus, sgn(A'P(Aj)) = {—iy^^sgn{Q{Xj)) where sgn(a) denotes the sign of a. Now if 
i-iy-^ sgn{Q{Xj)) = (-l)^(sgn(g(A,+i)), then - sgn(Q(A,)) = sgn(g(A,+i)), so Q{s) must 
change sign on the interval [Aj+i, A^]. Since Q{s) is a polynomial of degree r, and r < n, Q{s) 
can switch signs at most r times, else Q{s) = 0. But this means there are at least n — r — 1 
intervals [A^-^+i, A^J, for A; = 1, . . . , n - r - 1, for which sgn{Q{XjJ) = sgn{Q{Xj^+i)), and 
therefore sgn(KP(Ajj,)) = —sgp.{KP{Xj^^i)). So KP{s) must change sign over at least 
n — r — 1 intervals, and therefore has at least ri — r — 1 zeros on [A^, Ai]. Again, since the 
error is not identically zero when r < n, and the degree of KP{s) is n — r — 1, this implies 
that all the zeros of KP{s) lie in (— oo,0). 

Suppose with some p < r, Xi. = X^ for j = 1, . . . ,p. Observe from partial fraction expansions 
of II{s) and IIr{s) that the error can be written as a rational function of degree n + r — p—1 



over degree n + r —p with distinct poles, n of these poles belong to H{s) and the remaining 
r — p come from the poles of H^{s) that are distinct from the poles of H{s). Now let 

2r n— r— p— 1 n r—p 

R{s) = l[{s-s,), P{s)= n {s + z.^, Q(s) = l[is~K), g(s) = J](s-AfcJ, 

i=l j=l j=l /=1 

where {Afc,}[rf = {Xk}l.=i\{^i}7=i- Hence, H{s) — Hr{s) = K \^'~\^' . Observe that there are 

at most 2p subintervals of the form [Aj*, Aj.+i] or [Aj.-i, Aj.], where Aj. G {Xi}^=i fl {Afc}^=i. 
It follows that there are at least n — 2p — l subintervals of the form [Aj, Aj+i], where Aj, Aj+i ^ 
{Ai}"=i n {Afcj^^p On each such subinterval for which this is the case, we have 

Res{H{s) - Hris); A,) = k_JV^jM>^ = 0, > 0. (11) 

n (A. - A,)Q(A.) 

So sgn(ii'P(Aj)) = (— 1)*^^ sgn((5(Aj)). By the same argument as above where the poles of 
H{s) and Hr{s) are distinct, either Q{s) or P{s) has a zero on the interval [Aj, Aj+i]. Since 
Q{s) has at most r—p zeros, this means that there are at least n—2p—l — {r—p) = n—p—r — 1 
subintervals between poles of H{s) where P{s) has zeros. Hence, the lemma is proved. D 

Lemma 3.2. Let H{s) = b'^{sl — A)^^b be SSS, and Hr{s) = b^{slr — Ar)~^br be any 
reduced order model of H{s) constructed by a compression of H{s), i.e., Ar=QjArQr, 
br= Q^b. Then for any s>0, H{s) - Hr{s) > 0. 

Proof. Pick any s > 0. Then (s/„ — A) is symmetric, positive definite and has a Cholesky 
decomposition, {sin — A) = LL'^ . Define Zj. = L^Qr. Then 



H{S) - Hris) = b^ [{sin - A)-' - Qr {Ql{sln - A)Q,) Q 

= {L-'bY \l - Z, [ZjZ,)-' Zj] {L-'b). 



Note the last bracketed expression is an orthogonal projector onto Rar\{Zr)^, hence is positive 
semidefinite and the conclusion follows. D 

Our convergence analysis of IRKA will use its formulation as a fixed-point iteration. The 



analysis will build on the framework of [ll|. Let 



H{s) = y-!^ and Hr{s) = y^^ (12) 

be the partial fraction decompositions of H{s), and Hr{s), respectively. Given a set of r 
interpolation points {sj}^^^, identify the set with a vector s = [si, . . . , s^]^. Construct an 



interpolatory reduced order model Hr{s) from s as in Theorem 11.11 and identify {Xi}l^i with 
a vector A = [Ai, . . . , A^]^. Then define the function A : C^ — )■ C" by A(s) = —A. Aside 
from ordering issues, this function is well defined, and the IRKA iteration converges when 
A(s) = s. Thus convergence of IRKA is equivalent to convergence of a fixed point iteration 
on the function A(s). Similar to s and A, let = [0i, . . . , 0r]^- Having identified Hr{s) with 
its poles and residues, the optimal ^2 model reduction problem may be formulated in terms 
of minimizing the cost function J7'(</), A(s)) = \\H ~ -f^rll^j ' where 



J{4>, A(s)) = J2 HH{\) - Hr{Xi)) + 



i=l 






(13) 



See [Sl] for a derivation of ( TT3l) . Define the matrices Su, S12, S22 G 



[Su]u = -C\ + >^jr\ [s 



12 



«j 



-(A, + Aj)-2 and [S; 



22 



«j 



as 

-2(A, + A,; 



for i,j = l,...,r. Also, define R, E e MJ'^'': 

R = diag({0i, . . . , 0,}), and E = dmg{{H" {-~X,) - i/;(-Ai), . . . , H"{-Xr) - //^(-A,)}. 

Lemma 3.3. Let H{s) be SSS and let Hr{s) be an IRKA interpolant. Then E is positive 
definite at any fixed point of X{s). 



Proof. By Lemma |3.2[ H{s) — Hr{s) > for all s G [0,oo). Thus the points H(—Xi) — 
Hr{—Xi) are local minima of H{s) — Hr{s) on [0,oo) for z = 1, . . . ,r. It then follows that 
H"{—Xi) — H'^{—Xi) > 0. But by Lemma ISTT] H{s) — Hr{s) has exactly 2r zeros in C+, so 
H'i-Xi) - H';{-Xi) > for 2 = 1, . . . , r. D 



Lemma 3.4. The matrix S 



»-'ll >-'l2 
'S'12 >S'22 



is positive definite. 



Proof. We will show that for any non-zero vector z = [zi, Z2, . . . , Z2r\ G 



T r- n2r 



z'Sz 



E-^^''*-<E 



4 = 1 



■^■r+jS 



Ait 



1 2 



i=l 



dt > 0. 



Define Zr = [zi, Z2, ■ ■ ■ , Zr] G M.'^ and Z2r = [zr+i, Zr+2, • • • , Z2r] G M'". Then 



Z bZ — Z biiZj. + Z,Z Si2Z2r + -22j.'-'22'22 



(14) 



Let A = diag(Ai, . . . , A^) and it be a vector of r ones. Note that Su solves the Lyapunov 
equation ASu + SuA + uu^ = 0. Thus, 



z SiiZj. 



zJe^^uu^e^^Zr dt 



E 



1=1 



ZiC 



Xrt 



dt 



(15) 



Similarly, S'12 solves AS'12 + >S'i2A — Su = 0. An application of integration by parts gives: 



i=l 



i=l 



-100 / /"°° \ 



-Z S12Z2T 



(16) 



Finally, note that S22 solves AS22 + S'22 A — 2Si2 = 0. Repeated applications of integration 
by parts then yields the equality: 



Z2rS22Z2r= / t^ (J^ ^r+lC^') ^^ 



(17) 



j=l 



Combining equations dHD, ( TT5l) . ( TT6l) . and ( TT71) gives the desired results since 



^Sz = /; 



i=l ^ i=l 



dt. 



n 



Then it follows that the Schur complement ^22 — S12S11S12 of S is also positive definite. 
With the setup above, we may now commence with the proof of Theorem 13.11 



Proof of Theorein l3.lt It suffices to show that for any fixed point which is a local minimizer 
of J'{(p, A(V)), the eigenvalues of the Jacobian of A(s) are bounded in magnitude by 1. As 
shown in [ll|, the Jacobian of A(s) can be written as —S^^K where 



Sc = S22 — S12S11S12 and K = ER 



-1 



First off, note that from Lemma [3. 3 [ and the fact that H{s) is a ZIP system by Lemma [2. H 
K is positive definite. Evaluating the pencil K — XSc at A = 1 gives 



$ 



-022 + ER + iSi2'-'ii 012, 



This pencil is regular since Sc is positive definite by Lemma [231 and therefore det(l^ — XSc) 
is zero if and only if det(S'~^l^ — XI) = 0. 



Let V^tT" denote the Hessian of the cost function jT{4>,X{s)). As shown in ll|], V^JT" can 

be written as 



V'J 



I 





M 


I 


0" 


[0 


R\ 




[u 


r\ 



where M 



Sii S12 

i5i9 099 — ER 



Note that — $ is the Schur complement of M. Hence, if the fixed point is a local minimimum, 
then — $ must be positive definite and so for A = 1 the pencil is negative definite. Since 
both K and Sc are positive definite, there exists a nonsingular transformation Z by which 
the quadratic form y {K — XSc)y is transformed into z {A — XI)z, where A is a diagonal 
matrix formed from the solutions of 

det{K - XS,) = 0. (18) 

Thus, the solutions of fITS]) correspond to the eigenvalues of S~^K. A — I must be negative 
definite since $ is, and therefore the eigenvalues of the S~^K must be real-valued and less 
than one. Furthermore, note that P = S~^K solves the Lyapunov equation 

PS;' + s;'p^ = 2S;'KS;\ 

so by the standard inertia result, all the eigenvalues of S~^K are positive, and the desired 
result follows. D 
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